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ABSTRACT 

We numerically investigate the possibility of planetesimal accretion in circumbi- 
nary disks, under the coupled influence of both stars' secular perturbations and friction 
due to the gaseous component of the protoplanetary disk. We focus on one crucial pa- 
rameter: the distribution of encounter velocities between planetesimals in the 0.5 to 
100km size range. An extended range of binary systems with differing orbital param- 
eters is explored. The resulting encounter velocities are compared to the threshold 
velocities below which the net outcome of a collision is accumulation into a larger 
body instead of mass erosion. For each binary configuration, we derive the critical 
radial distance from the binary barycenter beyond which planetesimal accretion is 
possible. This critical radial distance is smallest for equal-mass binaries on almost cir- 
cular orbits. It shifts to larger values for increasing eccentricities and decreasing mass 
ratio. The importance of the planetesimals' orbital alignments of planetesimals due to 
gas drag effects is discussed. 
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1 INTRODUCTION 

At present more than 40 planets have been found in binary 
star systems, which corresponds to abo ut 20 % of the whole 
samp le of known extrasolar planets ijPesidera fc Barbieril 
l2007l ). All of them are on so— called S— type or circumprimary 
orbits that encircle one component of the double star sys- 
tem. Planets may also revolve on dynamically stable orbits 
about both binary stars on so-called P-type or circumbinary 
orbits (e.g. Holman and Wiegert 1999). At present, there is 
only one planet assumed to revolve on a circumbinary or- 
bit. The primary, HD 202206, is a metal rich main sequence 
star. Its companion, HD 202206b, is assumed to be a low- 
mass brown-dwarf with an estimated mas s of 17.5 Mjup ac- 
cording to radial ve locity measurements (|Udrv et al.ll2003 : 
ICorreia et al] l200d ). This mass is beyond the widely ac- 
cepted brown-dwarf limit of 10 Jupiter masses. The third 
body, HD 202206c, is assumed to be a Jupiter-like planet 
revolving about the binary. Adopting a mass of 2.41 Mjup 
for HD 202206c, Correia et al. (2005) show in a dynamical 
analysis that HD 202206c and HD 202206b might be in a 5/1 
mean motion resonance. This resonant configuration is pos- 
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sibly a consequence of an inward migration of HD 202206c 
due to the forces exerted by a viscous circumbinary disk 
(Nelson 2003). 

The signature of circumbinary planets in radial velocity 
measurements is in most cases difficult to be spotted because 
of the short-term large-amplitude velocity fiuctuations on 
the primary induced by the companion star. However, plan- 
ets might be very frequent among the population of close 
star couples. Mid-infrared emissions have in fact revealed 
substantial circumbinary material around PMS close bina- 
ries like DQ Tau, UZ Tan, GW Or i (Mathieu et al. 2000) and 
AK Sco (|Jensen fc Mathieulll99'it). and more detached bina- 
ries l ike GG Tau l|Dutrev et al■lll994^ ■ UY Aur (|Close et al.l 
|l998l ). These disks are significantly more massive than the 
minimum-mass solar nebula which suggests abundant plan- 
etary formation. However, both concurrent mechanisms for 
extrasolar planet formation, namely the standard solid- 
core model forming planetary embryos fr om runaway and 
oliga r chic accretion of planetes i mals (e.g iGreenberg et al.l 
Il978l : IWetheriU fc Stewartlll989l : iKokubo fc Idall200ol) pos^ 
sibly followed by gas infall onto the core for the formation 
of giant planets (Pollack et al. 1996; Bodenheimer and Lin 
2002), and the alternative model of local gravitational col- 
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lapse of disk material (|Bosslll997l ) might be affected by the 
secular perturbations of the binary. 

We will focus here specifically on the "standard" solid- 
core scenario, in particular on the stage of planetary em- 
bryo accretion. This stage is just before the one investigated 
by Quintana and Lissauer (2006) who determined regions 
around close binaries where planets can form by accumu- 
lation of embryos. Within the frame of this model, the ini- 
tial stage of planetesimal accretion is the runaway growth 
of isolated planetary embryos. This stage is very fast and 
efficient, provided that the encounter velocities among plan- 
etesimal are low, i.e. muc h smaller than the escape velocities 
of the growing embryos (|Lissaueij|l993) . As a consequence, 
runaway growth is particularly sensitive to any external per- 
turbation able to stir up relative velocities in the planetesi- 
mal swarm. In the specific case of a circumbinary disk, rel- 
ative velocities between planetesimals can be expected to 
be higher as compared to the single star case. In the vicin- 
ity of the binary, velocities could be increased to values for 
which runaway accretion is no longer possible, or even to 
higher values for which collisions may result in mass removal 
via cratering or fragmentation rather than accretion. The 
threshold velocity for the stop of runaway is of the order 
of the escape velocities iiesc(ij) of the growing bodies while 
the threshold velocity Vero for the limit between accreting 
and eroding impacts is less straightforward to estimate and 
depends on the respective sizes and physical properties of 
the two colliding planetesimals. In any case, the distribu- 
tion of encounter velocities is the key parameter controlling 
the fate of a coUisionally interacting swarm of planetesimals. 
This is the issue which we address here, by numerically esti- 
mating the distribution of < At; > in a disk of planetesimals 
surrounding a binary system in order to derive assertive con- 
clusions on the possibility of planet formation. 

Only taking into account the gravitational potential of 
both stars might lead to incorrect estimates of relative ve- 
locities, in particular close to the star couple. Indeed, the 
planetesimal accretion phase is believed to take place while 
a large fraction of the primordial gas disk is still pr esent . Gas 
drag is known to circularize orbits ( Adachi et al. 19761 ) and 
thus reduce relative vel ocities in dynamically hot systems 
IIMarzari fc ScholllbOOd ). Besides circularization, gas drag 
has the additional effect of aligning orbital ap sides, an effect 
which is also called orb ital synchronization l|Marzari et al.l 
Il997l : iMarzari fc ScholfcoOO l . The consequences of this syn- 
chronization on relative velocities depend on planetesimal 
sizes. Bodies of the same size all have their periapses aligned 
in the same direction. Smaller bodies have a much smaller 
dispersion of periapse directions than larger bodies. Further- 
more, the mean direction of the peria stron alignment is als o 
different for objects of different sizes (jThebault et al.ll2006h . 
As a result, while relative velocities between small planetes- 
imals of the same size may drop to almost zero, relative 
velocities among bigger bodies are more significant because 
of a larger periapse direction dispersion. The highest relative 
velocities can be expected between objects of different sizes 
due to the size-dependency of the periastron alignment. The 
relative velocity distribution depends, therefore, not only on 
the radial distance to the barycenter of the binary but also 
on the sizes of colliding planetesimals and on the gas density. 

In a pioneering attempt to model planetary accretion in 
circumbinary disks, Moriwaki and Nakagawa (2004) used a 



purely gravitational model without gas drag. They found 
very stringent limits on the minimum distance from the 
binary beyond which planetesimal accumulation is possi- 
ble. However, as pointed out above, neglecting the velocity 
damping effect and orbit alignment due to gas friction may 
yield misleading results for the efficiency of planetesimal ac- 
cumulation, in particular in the vicinity of the binary. In the 
present work, we redress this issue by performing numerical 
simulations where the effects of gas drag are taken into ac- 
count, although with unavoidable simplifying assumptions. 
Moreover, we compute the average encounter velocities be- 
tween planetesimals from their tra jectories and not from 
the va lues of forced eccentricities as iMoriwaki fc Nakagawal 
(j2004i ) who used the relation Aw ~ e.Vkep, where Vkep de- 
notes Keplerian velocity. This often used relation holds only 
under the condition that all orbits have fully randomized 
Keplerian angles. However, in circumbinary disks with gas, 
orbital phasing of periastra directions is of high importance 
because of the combined effects of secular perturbations and 
gas friction. Another limitation of the Av ~ e.Vkep relation 
is that it works locally and it does not take into account 
the possibility of large radial excursions of bodies with high 
eccentricitie s (for more on this subject, see for instance the 
discussion in lThebault fc Doressoundiramll2003l ). Both these 
issues suggest that the relation used by Moriwaki and Nak- 
agawa (2004) may not be appropriate in presence of large 
orbital eccentricities and possibly of gas drag. Our numerical 
approach, which takes into account all close particle encoun- 
ters representing collisions, has the advantage to handle au- 
tomatically orbital phasing (and dephasing) as well as radial 
mixing effects. This allows us to derive more reliable values 
for the limiting distance from the barycenter of the binary 
where planet formation is possible for different parameters 
of the binary. Not surprisingly, our results significantly differ 
from those of Moriwaki and Nakagawa (2004) in particular 
close to the binary where planetary accretion is more likely 
to form planets due to the higher density of material. 

Due to obvious numerical constraints, mutual gravita- 
tional interactions among planetesimals, like for instance 
dynamical friction, are not included even if, in principle, 
they might be taken into account in N-Body simulations. 
However, the speed of current computers is still far too low 
to allow a systematic investigation of binary systems sur- 
rounded by a self-gravitating planetesimal disk. These sim- 
ulations are still restricted to disks with massless planetesi- 
mals. Their size distribution, of course, does not evolve with 
time. To do so, particle-in-a-box computations, which can 
follow the evolution of a large number of massive bodies, 
are required, but the price to pay is a significant loss of ac- 
curacy in the computations of the orbital evolution of the 
planetesimals that, in the present scenario, is crucial. As 
a consequence, we decided to adopt a deterministic numeri- 
cal approach, which, although neglecting important physical 
effects, gives a reliable estimate of encounter velocity distri- 
butions and thus enable to derive "maps" of the disk regions 
with relative velocities low enough to favor planetesimal ac- 
cretion and planet formation. 

Our numerical model is described in detail in Sec- 
tion 2. Section 3 is devoted to the computations per- 
formed in the gas-free case, a scenario similar to that of 
iMoriwaki fc Nakagawal (|2004l ). In Section 4 we use the re- 
sults obtained from models including gas drag in order to 
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map the regions where we expect planet formation by accre- 
tion. Section 5 is dedicated to the dicussion of our results. 



2 THE NUMERICAL MODEL 
2.1 Main characterstics 

Our deterministic code is an upgraded version from the on e 
used in several earlier studies l|Marzari fc Scholilll998l . [20001 ) . 
adapted to the specificities of the circumbinary configura- 
tion. We consider a system of A*' ~ 10000 planetesimals, 
treated as non-gravitationally interacting test particles, sub- 
mitted to the gravitatio nal potentials of both stars and to a 
gas drag force following lAdachi et al.l (|l976l ) and described 
in the following section. Planetesimal orbits are integrated 
using a fourth-order Runge-Kutta scheme with fixed step 
size. We like to point out that for the circumprimary prob- 
lem (|Thebault et al.ll2006l ) we applied a different numerical 
scheme which is less convenient for the circumbinary case. 

The binary stars move on an elliptic orbit of eccentricity 
Bb and semi-major axis at- For sake of simplicity and to limit 
the number of free parameters to a manageable value, we as- 
sume for all runs ai, = 1 AU. This is the standard value also 
taken by iMoriwaki fc Nakagawal l|2004l ) in their numerical 
explorations. We are thus here implicitly focusing on plan- 
etesimal accretion in the regions typically beyond 3-4 AU 
from the binary's centre of mass. This region is sufficiently 
away from the critical instability region found by Holman 
and Wiegert (1999) to guarantee orbital stability for almost 
planar and circular motion. Both stars' total mass is fixed 
and equal to 1 solar mass. The primary parameters of our 
model are thus the mass ratio and eccentricity of the binary 
star system. We use five different mass ratios q for the bi- 
nary, q = 0.5, 0.4, 0.3, 0.2 and 0.1. The stellar masses mi of 
the primary and m2 of the companion are then mi — 1 — q 
and m2 — q, in solar mass units. For each mass ratio, models 
with six different binary eccentricities ei, are simulated with 
Bb = 0.0, 0.1, 0.2, 0.3, 0.4 and 0.5. This gives a total of 30 
models. Each model is simulated with and without gas drag. 
In models taking into account gas drag, we consider plan- 
etesimal of physical radii ranging from 0.5km to 100km. The 
initial eccentricities and orbital inclinations are chosen with 
a uniform random distribution between and 10~^. The 
initial particle semi-major axis are randomly distributed 
between amin = 4AU and amax ~ 12 AU. The value for 
dmin being approximately equal to the inner truncation of 
a circumbinary disk by the tidal influence of the binary 
l|Artvmowicz et al.|[l993 : ICunther fc Klevll2002l '). Whenever 
the semimajor axis of a planetesimal decreases below the 
inner limit because of drag friction, the body is removed. 

The typical timescale for one run is ~100 000 binary 
revolutions, which is approximately the time necessary to 
reach a steady state, with fully developed secular perturba- 
tions and the gas drag induced orbital phasing. After such 
a steady state is reached, we determine relative velocities 
among the planetesimal population. For this purpose, we use 
the positions and velocities of all planetesimals at a fixed 
time. For each planetesimal, we search its nearest neigh- 
bouring planetesimal and compute the relative velocity cor- 
rected for Kepler shear due the different radial distances. 
Relative velocities are sampled in bins in radial direction. A 



bin size of 0.5 AU is used. The median value for each bin 
is then taken as the representative value for the coUisional 
or impact velocity of planetesimals at the corresponding ra- 
dial distance. Since, as pointed out, the orbit alignment is 
size dependent, we compute the median relative velocities 
^'V{Ri,R2) between planetesimal populations for each pairs 
of radii Ri, R2 in each bin. We typically use 20 to 30 close 
encounters in a bin to determine the corresponding median 
relative velocity. 



2.2 The gas drag force 

The drag force is modelled in laminar gas approximation by 

= -iCUreiVrel, (1) 

where K is the drag parameter given by (|Karv et al.lll993h 



A- = 3 



PgCd 
8ppis 



(2) 



where s is the planetesimal radius, ppi its mass density, pg 
the gas density of the protoplanetary disk and Cd a dimen- 
sionless drag coefficient related to the shape of the body 
(~ 0.4 for spherical bodies). Exploring the gas drag density 
profile as a free parameter would be too CPU time con- 
suming and we shall restrict ourselves to one pg distribu- 
tion. We assu me the st a ndard Minimum Mass Solar Neb- 
ula (MMSN) of lHavashil (|l98ll ). with pg = pgo(a/lAU)~'^-'^^ 
and Pgo = 1.4 x 10~^g.cm~^. We take a typical value 
Ppi = 3g.cm~^. 

We are here implicitly assuming that the gas disk is ax- 
isymmetric and pressure supported. We are aware that this 
axisymmetric prescription does probably not represent the 
real physical structure of the circumbinary gaseous compo- 
nent of the disk. This is a complex issue, which has only 
been addres sed in a handful of elaborated numerical investi- 
gations fe.g lArtvmowicz et al. lll994l : lGunther fc Klevl[200^ : 
lNelsonll2003l : |Pichardo et al.ll2005l ). All these studies clearly 
show the inner tidal truncation of the disk, but are less de- 
tailed on the structure within the disk, which is particularly 
complex to be precisely computed, especially in the eccen- 
tric binary case. The only clear result is the onset of complex 
time-evolving spiral patterns extending from the inner edge 
of the disk towards the outer regions. At present it appears 
a very demanding task to perform a fully coupled numerical 
treatment of both the gas disk and an embedded planetesi- 
mal swarm. For this reason, we prefer in a first step a simpli- 
fied approach where the gas drag force is just described by 
Equ[T] The onset of density waves in the disk induces radial 
and azimuthal variations both of the disk density and of the 
gas velocity where an additional radial component adds up 
to the Keplerian tangential component. These effects might 
alter the strength of the drag acceleration when a planetes- 
imal enters and exits the higher or lower density region. 
However, it is reasonable to assume that, on the average, 
this effect would not strongly modify accelerations obtained 
by Equ[T]and thus the dynamical evolution of planetesimals, 
at least for planetesimals in the km-size range. In particular, 
we expect the occurrence of the periastron alignment. 

As already mentioned, we start all planetesimals on al- 
most circular orbits. This choice is from a dynamical point 
a view equivalent to assuming that, at t = 0, the proper 
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eccentricity is almost equal to the forced one ef. This is 
a somewhat artificial and not a necessarily realistic start- 
ing configuration motivated by the fact that, at present, we 
have no generally accepted clues either on how planetesimals 
form within a disk with spiral waves nor on the initial or- 
bital elements of the planetesimals as they detach from the 
gas. We integrate the model until a steady state is reached 
where periapses are aligned and the interaction with the gas 
becomes only a minor perturbation of the Keplerian orbit. 
This procedure is valid for small planetesimals, since they 
reach the same periapse alignment independently of their 
starting values with the same distribution of proper eccen- 
trcities. Larger planetesimals, which are much less affected 
by gas drag, have their proper eccentricities increasing dur- 
ing their growth by collisions and close encounters with other 
large planetesimals. As a consequence, they too should lose 
memory of their initial orbits. As a consequence, we think 
that our initial choice of eccentricities does not introduce 
too significant effects on the calculation of relative impact 
velocities once the steady state is reached. 



Analytical fit - 
N-body simulations 
MorlwakI and Nakagawa (2004) - 
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Figure 1. Relative impact velocities between planetesimals com- 
puted with ou r equation [3] (continuous lin e) compared to the 
predictions of llMoriwaki fc Nakagawal |2004) (dashed line). The 
empty circles show the numerical results of the N-body integra- 
tions. The parameters of the binary are q = 0.2 and ej, = 0.4. 



3 IMPACT VELOCITIES IN THE GAS FREE 
CASE 

In the gas free models we use 25000 particles with semima- 
jor axes ranging from 3 to 50 AU. However, we focus here 
essentially on the region within 12AU, since, as it clearly 
appears in Fig[T] average encounter velocities drop down to 
almost zero beyond this limit. For the r < 12 AU region, 
the (Av) we obtain might be fitted, for the range of stellar 
masses are orbital eccentricities explored, by the empirical 
formula 

{Aw)(^) = (-384.2 -h 2067.7a - 447522efl 
-3864.9q^ -h9462.4gefl) 
-(-(1024.8 - 5783.4 + 16431es 

4-11251?^ - 34396qei3)-^ (3) 

which is a second-order polynomial in q and et multiplied by 
an additional term equal to y/{r) accounting for the radial 
dependence of (Av) on the local Keplerian velocity. 

Since we neglect the gas drag force, the radius of the 
planetesimals is not a parameter of the problem. In Fig[T] 
we also com pare our resu l ts to the (Av) computed with 
theiMoriwak i fc Nakagawal (|2004l ) analytical approximation 
for the specific case with q = 0.2 and e;, = 0.4. For ref- 
erence, we also plot the numerical data obtained from the 
N-body integration. It can be clearly seen that our rmmer- 
ically derived values are always significantly smaller than 
the|Moriwaki & Nakagawa (2004) estimates. This shows the 
weakness of the approximate relation An ~ evhe.p for the 
present case, where orbital phasing and radial excursions of 
the planetesimals play a crucial role in determining the av- 
erage impact velocity. As an example, while our simulations 
give low relative velocities, of the order of 50 m/s at about 
11 AU from the baryce n ter, th e analytical approximation of 
iMoriwaki fc Nakagawal ()2004h predicts such low velocities 
only beyond 50 AU. Our numerical results give much lower 
values for impact velocities which allow planetary formation 



to occur much closer to the barycenter of the binary even in 
absence of gas friction. This is a first indication that plan- 
ets on circumbinary orbits may be abundant due to the low 
relative velocities in the regions near the binary with higher 
surface density. In the following section, we show how this 
promising result is affected by gas drag. 



4 IMPACT VELOCITIES IN THE PRESENCE 
OF GAS 

We now turn to the more realistic case including gas drag. 
As previously mentioned, it is well known that the combina- 
tion of the drag force and perturbations by the gravity field 
of the binary causes a strong paeriapse phasing of the orbits 
of the planetesimals. The planetesimals' physical sizes are 
now a crucial parameter. Bodies of similar sizes have their 
periapses aligned towards approximately the same direction, 
which leads to very low impact velocities. Different size plan- 
etesimals are aligned towards different directions, increasing 
the non-tangential component of the impact velocities. 

To compute the relative velocity between any possible 
pair of different sized planetesimals in the swarm is far be- 
yond the present computing capabilities. Therefore, we de- 
rive the mean relative velocities Ad for two separate cases 
assumed to be representative for the accreting capability of 
a swarm: 

• A«jji,_Ri between equal-size bodies of radius K\. In this 
case the periapse alignment is stronger. 

• Avfi\,Ri for a target with radius R\ and a projectile 
with radius Ri — Ri/2. This size ratio has been chosen since 
it is the configuration for which, in the la rge majority of 
cases explored in our circumprimary study (iThebault et al.l 
[200tf ). the kinetic energy delivered by the impactor is maxi- 
mal. For Ri/2 < R2 < Ri the relative velocity Avri^r^ de- 
creases and approaches the low velocity value of equal-size 
impacting bodies. For R2 < Ri/2, the increase in Aw_Rj,_r2 is 
too slow to compensate for the diminishing impactor mass. 

For each of the 30 different {q,et) binary configura- 
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Figure 2. Figure showing the relative velocity between planetesimals for different values of the binary eccentricity ej,, mass ratio q and 
planetesimal sizes R\ and R2. The plots are organized as follows: top left q = 0.5 (equal-mass stars case), Hi = 25 km and R2 = 50km; 
top right q = 0.1, _Ri = 25 km and R2 = 50km; bottom left q = 0.1, Ri = 0.25 km and R2 = 0.5km; bottom right q = 0.1, Ri = 100 km 
and -R2 = 50km. 



tions considered, we perform 6 simulations exploring differ- 
ent target-impactor pairs (J?i,l/2i?i) in the 0.5km< _Ri < 
100km range (leading to a total of 180 runs). We compute 
for each pair the average impact velocity < Av >ri^r2- In 
Fig[2] we show some examples of the relative velocity be- 
tween different size planetesimal swarms. A general feature 
is that binaries with higher eccentricities always lead to 
higher {Av)m^R2- Conversely, the most perturbing config- 
urations are the ones with lowest q values, the equal-mass 
case {q = 0.5) leading to the lowest {Av)ri^r2- This result is 
logical when considering the structure of the secular term of 
the perturbing potential (see for example the Appendix of 
iMoriwaki fc Nakagawa 2004): the forced eccentricity term is 
proportional to (1 — 2g), i.e. it vanishes for the equal-mass 
case (where only short period perturbation terms remain). 
An alternative and simplified, but partially correct explana- 
tion is that the equal mass case is the one for which the radial 
excursion of both stars from the centre of mass is the more 
limited, i.e., they both stay at 0.5 AU from it, whereas for 
smaller 7712 /mi values, star number 2 moves further away 
from the centre of mass and "enters" deeper into the cir- 
cumbinary disk. Within the range of parameters which have 
been explored, the most perturbed case is thus logically the 
one with mass ratio g = 0.1 and et = 0.5. In any case, even 
for different size impacting objects, one robust result is that 



the presence of gas drag always strongly reduces the relative 
velocity between planetesimals. This can be clearly seen in 
Table I, where we display, for one given set of binary param- 
eters (the same as in Fig[T]), relative velocity values for the 
case without gas drag, both from our simulations and from 
Moriwaki and Nakagawa (2004), and for the case with gas 
drag for a target planetesimal of size _Ri = 5km. We consider 
both the case of equal-size planetesimals with R2 = Ri = 5 
km and different size planetesimals with Ri — 5km and 
R2 = 2.5 km. 

As discussed in the introduction, all {Av)ri^r2 values 
have to be compared to two critical values in order to de- 
termine the coUisional evolution of a planetesimal swarm. 
One is the escape velocity of the target v^gc(Ri) , which can 
easily be derived from the body's radius and density. Only 
if {Av)ri^r2 « Uesc{Hi) can runaway growth occur. The 
other crucial threshold value is Vero, the threshold velocity 
beyond which the net outcome of an impact is mass ero- 
sion instead of accumulation into a larger body. To estimate 
the value of «ero{i{i,fl2) for a given projectile and target is 
a complex task. It has to be derived from models of cra- 
tering and fragmentation which, at present, give uncertain 
results due to our poor k nowledge of the si ze-st rength scal- 
ing law. As discussed in iThebault et al] (|2006l ). there are 
many different models in the literature giving the value of 
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the critical specific shattering energy Q* and the size distri- 
bution of fragments as a function of the projectile and target 
size, of composition a nd of impact vel o city. T o be conserva- 
tive, we consider as in iThebault eTsU l|2006l) thr ee different 
estima t es of Vero(Ri. r 2 \ corn pute d according to Holsappl3 
j 19941 ). iMarzari et ^ (119951') and ISenz fc Asphaud (|l999l f 
prescriptions for Q* and IThebault et al. I (|2003l) prescription 
for the fragmented and cratered masses. Out of these we 
take the maximum value that gives an optimistic estimate 
of the region where accretion is possible, and the minimum 
and more conservative value. 

When we compare {Av)ri,r2 to v^sc{ri) and «ero(fli,fl2) 
we can outline 3 different accretional modes: 

• 1) (Av) « iiesc(fli)- In this case, runaway accretion 
can proceed as in the 'standard' unperturbed scenario. The 
additional perturbations of the binary are not large enough 
to cancel out the gravitational focusing factor {vesc/ {Av))'^ 
of t he target, whic h is the source of this fast growth mode 
e.g. (|Lissaueij|l993l ). 

• 2) For values of (Av) larger than v^sc but still smaller 
than Vera, planetary accretion is still possible but runaway 
growth will either not occur at all or it will start only when 
large bodies have formed such that Vesc > (Av) . In this case 
planetesimals will accumulate into planetary bodies but at 
a significantly slower rate. 

• 3) If (Av) exceeds Vero, the majority of collisions end up 
with cratering and fragmentation that overcome accretion. 
The planetesimal population, instead of growing in size, will 
be slowly ground down to dust. 

For sake of clarity we display our numerical predictions 
on planetesimal accretion in 4 summarizing graphs. In each 
of them we display, for all explored values of q and et, the 
radial distance ri beyond which planetesimal accretion is 
possible, for a population of objects having an initial size 
-Ri = Rmin, all the way up to the maximal size of 100km 
considered in our simulations. Beyond that size, accretion 
dominates over erosion for all explored cases. In other words, 
n is the limit beyond which {Av)ri,r2 < Uero(Hi,H2) for all 
{Ri,R2) pairs in the _Ri > Rmin range. We consider two 
cases for the starting Rmin size: one "standard" case with 
Rmin = 5km, and one "small planetesimals" case for which 
Rmin ~ 0.5km. Note that the value Rmin = 5 km is the 
commonly accepted minimum size for planetesimals: Inde- 
pendently of their formation process, when they reach this 
size they detach from the gas of the disk and move on inde- 
pendent Keplerian orbits, only perturbed by the gas. This 
value is however approximate and this is why we also explore 
a lower value of 0.5km for the initial planetesimal size. We 
also consider the two maximum and minimum values derived 
from the different prescriptions considered for «ero{Bi.fl2)- 

Fig|3] displays the most favorable case for accretion, 
where we assume Rmin ~ 5km and adopt the highest, and 
thus less restrictive value for Vero- It can be seen that, in 
spite of the binary gravitational perturbations and of the 
periastra de-phasing of different size planetesimals, for mass 
ratios q —0.5 and 0.4 planet formation proceeds in all the 
regions of the disk and for all binary eccentricities. The only 
region where accretion is inhibited is within the inner tidal 
gap of the disk. For smaller values of q and high binary ec- 
centricities eb, the inner border for accretion shifts to larger 
radial distances. For q — 0.1 and et — 0.5 planet formation 




mass of the secondary 

Figure 3. Map of the radial distance r; beyond which planetary 
formation is possible, i.e. where all (i?i,i?2) encounters lead to 
mass accretion {(Ai')i{i i{2 < ^ero{ijl,ij2)) for ^-H objects of size 
Ri > Rmin = 5km, as a function of the binary mass ratio q and 
binary eccentricity ej. The color coding used in Figures 3-6 is the 
following: 

-green: ri ^ 4AU (the inner edge of our planetesimal disc) 
-pale blue: 4AU< n < 6AU 
-dark blue: 6AU^ r; ^ 9AU 
-red: 9AU< n < 12AU 

-black: ri ^ 12AU (the outer edge of our planetesimal disc) 
The radial distance given at the center of some rectangles is the 
minimum value beyond which runaway accretion is possible (i.e. 
(Ai')fli,H2 << Vesc(Ri)) for all _Ri > 5km. 



can occur only beyond 10 AU: strong secular perturbations 
due to the large eccentricity and low mass ratio of the binary 
prevent planetesimal accumulation closer to the barycenter. 
Note however that for most cases where we found accre- 
tion to be possible the growth mode can probably not be 
runaway, or at least not 100% runaway. This is because 
{Av)ri^r2 values are not << Vesc{fli) for all explored -Ri 
values. The only exceptions, where unperturbed runaway 
accretion is possible, are obtained for the equal-mass case 
q — 0.5 and low binary eccentricities ej, ^ 0.3 (the white 
rectangles in Fig [3}. 

For a more conservative choice of the erosion velocity 
(Fig|4]), the inner limit for planetesimal accretion moves sig- 
nificantly outward, as expected: only for q = 0.5, i.e. equal 
mass stars, is accretion possible at any radial distance, in- 
dependently of Bb. 

The cases with smaller 0.5km initial planetesimals ap- 
pear to be more critical ([S]and Fig[SJ. While for mass ratios 
of 0.5 and 0.4 planet formation is (as in the Rmin ~ 5km 
case) possible almost all over the disk, when q becomes 
smaller than 0.4 in most cases accretion of small planetes- 
imals is inhibited within 12 AU from the stars. The main 
reason for this accretion inhibiting behaviour is that for 
very small bodies (in the sub-kilometer range) , orbital phas- 
ing due to gas drag is very efficient and, as a consequence, 
the differential phasing between objects of different sizes is 
very strong. This leads to high (Aii)_ri,_r2, especially when 
compared to the very low erosion threshold velocities corre- 
sponding to such small objects. 
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(our simulations) 


(Moriwaki & Nakagawa 2004) 
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4.25 
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7.75 
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754m.s~l 


3 m.s^^ 


9m.s-l 


11.75 


49m.s~l 


444 m.s~^ 


5 m.s~^ 
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Table 1. Examples of impact velocity values < At) > obtained with different models, for the same example binary configuration as in 
Figll](5 = 0.2, ei, = 0.4). 




0,5 0.4 0.3 0.2 0.1 0.5 0.4 0.3 0.2 0.1 



mHss of the accoridary mass of the secondary 

Figure 4. Same as in Fig|3]but assuming the minimum, and thus Figure 5. Same as in Fig[3] ("optimistic" high Vero value) but 

more restrictive value of ti£^(,(fli^2) ■ assuming a minimum value of Rl equal to 0.5 km (small plan- 

etesimals). 



It is important to note that our results can be rescaled 
to different values of the reference gas density Pgo. Indeed, 
the behaviour of a population of planetesimals of size R 
in a gas disk of density pg is similar to the behaviour of 
a population of planetesimals of size R/X in a gas disk of 
density p^/X (see Equs.l&2). As an example, reducing by 
a factor 10 the size of the bodies is equivalent to model 
the evolution of larger bodies but increase the gas density 
by a factor 10. As a consequence, another way of reading 
Figs[S]and[S]is that in denser circumbinary disks, planetary 
formation may be easily inhibited even when starting from 
"standard" Rmin = 5km initial planetesimals. Note however 
that the rescaling is not fully straightforward, because 5km 
and 0.5km objects do not have the same threshold velocity 

Vero . 



5 CONCLUSIONS AND PERSPECTIVES 

With a full numerical approach we have investigated plan- 
etesimal accretion and the possibility of planet formation 
in circumbinary disks. The trajectories of 10000 planetesi- 
mals have been integrated over time. The evolution of one 
crucial parameter, i.e. the distribution of mutual encounter 
velocities within the system, is numerically explored for a 
wide range of binary parameters, specifically the mass ra- 
tio q and orbital eccentricity Eb, the semi-major axis being 
fixed at lAU. From these data we map the minimum radial 
distance from the barycenter of the two stars beyond which 
accretion is always possible as a function of the binary mass 



mass of the secondary 

Figure 6. Same as in Fig|4] (restrictive low Vero value), but as- 
suming a minimum value of Rl equal to 0.5 km (small planetesi- 
mals) . 

ratio and eccentricity. Planets can form only in regions of 
the disk beyond this distance. 

Our results significantly differ fr om those of previous 
studies (jMoriwaki fc Nakagawal |2004| ) mainly for two rea- 
sons: 

• The use of the dv oc e x VKep formula that relates the 
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forced eccentricity of the binary companion to the average 
collisional velocity does not work properly when there is a 
significant radial mixing of planetesimals. The comparison 
of our full numerical results with those obtained from the 
forced eccentricity differ by almost a factor two. The esti- 
mate based on the forced eccentricity leads to artificially 
higher impact velocities. 

• We include in our simulations the crucial effect of the 
drag force due to the gaseous component of the disk. This 
significantly affects the relative velocities between planetes- 
imals: the (At;) between similar-size bodies are significantly 
reduced by the alignment of the periapses, but different-size 
planetesimals are aligned towards different directions and 
they may experience higher impact velocities. 

Summarizing all our results on equal and different-size 
planetesimals we outline the regions of the disk where ac- 
cretion is always possible. Wc find that for a "standard" 
system where planetesimals have an initial size of~ 5km, 
accretion is always possible, in the considered > 4AU re- 
gion, for equal-mass binaries. For lower m2/mi mass rar 
tios, results significantly depend on the threshold velocity 
'Vero{Ri,R2) abovc wliicli ail impact by an impactor R2 on 
a target Ri leads to net erosion instead of net accretion. 
A general result is that lower m2/mi and higher ei values 
result in higher (At;)(/{i ^{2) and are thus less favorable to 
planetesimal accretion. 

One aspect of the problem that needs further investiga- 
tion is the initial orbital distribution of planetesimals. This 
parameter is linked to the difficult issue of when and how 
they detach from the gas. While fully formed planetesimals 
may be marginally affected by the spiral waves of the disk 
caused by the binary, smaller bodies would tend to follow 
the gas streamlines when they are small enough to be cou- 
pled to the gas. This is a complex problem since it is strongly 
tied to the planetesimal formation process. How is the dust 
aggregation process influenced by the gas waves and by the 
stellar secular perturbations? Is the gravitational instabil- 
ity still possible? These questions may be answered within 
a model capable to handle gas and dust at the same time. 
Such a complete model not being available yet, wc chose to 
make in our simulations the simplifying assumption that at 
t = planetesimal orbits are circulax, i.e., that their proper 
eccentricity is equal to the forced one. However, we think 
that this choice may not have too dramatic consequences 
on the estimation of average relative velocities. Indeed, we 
observe that the further evolution of the system, under the 
secular perturbations of the stars and of gas drag , leads to a 
steady state where the proper eccentricity is damped down 
and the periapses are aligned depending on body size. This 
steady state only weakly depends on the initial orbital pa- 
rameters, especially for smaller planetesimals. In addition, 
since the major contribution to impact velocities is due to 
the different alignment of different-size planetesimals, rather 
than due to the absolute value of eccentricity, another choice 
of initial planetesimal orbits may not lead to significantly 
different results. 

We considered a binary semimajor axis ai, of 1 AU in 
our simulations, but equation 3 for the gas-free case can be 
easily scaled to larger (or smaller) values of a^. However, 
when we include gas drag, such an easy scaling is no longer 
possible. The dynamical evolution of planetesimals depends 



on the local gas density, which in turns depends on the ra- 
dial distance to the star. As a consequence, the evolution of 
a planetesimal Ri at a distance r of the centre of mass of 
a binary of separation at is not equivalent to the evolution 
of the same planetesimal at distance r/2 around a binary of 
separation ab/2 (even when corrected by a timescale correct- 
ing factor), because the gas density is not the same at these 
2 locations. Furthermore, there is no reason to assume that 
the global gas disk profile is the same around binaries with 
different orbital separations. As an example, if for ab ~ 1 
AU it may be reasonable to adopt the minimum mass solar 
nebula density po at 4 AU, in a disk around a binary with 
semimajor axis of at = 5 AU, adopting the same density 
po at 20 AU appears inconsistent. In other words, while the 
relative velocity in absence of gas scales with at because of 
its dependence on the Keplerian velocity, its density does 
not necessarily scale with r when gas drag is considered. As 
a consequence, our results may be applied to disks around 
binaries with small separations, but they cannot be scaled 
to large at- 

An additional effect that we neglect in our simulations 
is the physical outcome of mutual collisions (be it accretion 
or erosion). As noted in the introduction, a model capable 
to handle at the same time the collisional and dynamical 
evolution of a swarm of planetesimals is still far from the 
present computing capabilities. 

Our work must be considered a step forward of a bet- 
ter understanding of planet formation in circumbinary disks. 
There are certainly aspects of the problems that need fur- 
ther investigation in future papers like planetary accretion 
in habitable regions axound binaries, planetary migration 
driven by disk interaction and by planetesimal encounters. 
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